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Abstract 

We derive minimal discrete models of the Boltzmann equation consistent with equilibrium ther- 
modynamics, and which recover correct hydrodynamics in arbitrary dimensions. A simple analyti- 
cal procedure of constructing the equilibrium for the nonisothermal hydrodynamics is established. 
A new discrete velocity model is proposed for the simulation of the Navier-Stokes-Fourier equation 
and is tested in the set up of Taylor vortex flow. For the lattice Boltzmann method of isothermal 
hydrodynamics, the explicit analytical form of the equilibrium distribution is presented. 
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Minimal kinetic models, and primarily the lattice Boltzmann method (hereafter LBM), 
have recently met with significant success in simulations of complex isothermal hydrody- 
namic phenomena. A few examples of the successful application are the simulation of fluid- 
particle suspensions, turbulent flows, spinodal decomposition M, ^ The first large-scale 
simulations of 3D spinodal decomposition in inertial regime and simulation of Brown- 
ian short-time regime j|] are a few successes achieved in the field of computational fluid 
mechanics through use of these approaches. The ability to handle a very complicated ge- 
ometry in a very simple manner has allowed these method to emerge as an alternative to 
purely continuum approaches, even for solving isothermal Navier-Stokes equation ||. In 
these methods, hydrodynamic equations (for example, the isothermal Navier-Stokes and the 
nonisothermal Navier-Stokes-Fourier equations in our terminology) are not addressed by a 
direct discretization procedure. Instead a simplified kinetic equation is introduced in such 
a way that hydrodynamic equations are obtained as its large-scale long-time limit. Two 
central issues in the construction of such models are the choice of discrete velocities, as few 
as possible, and the construction of the local equilibria such that the desired hydrodynamic 
equations are reproduced as closely as possible by the kinetic model. 

Unlike the isothermal case, the kinetic modeling of the nonisothermal hydrodynamics is 
a hitherto unsolved problem [|T], The kinetic models with proper thermodynamics are 
especially needed for the simulation of chemically reactive flows and the multiphase flows 
and near continuum flows in microdevices, which are difficult to simulate using a purely 
continuum models. However, the nonisothermal model is not established even in the case of 
the Navier-Stokes-Fourier dynamics of a single fluid. Apart from the conserved moments of 
distribution function, the mass, the momentum and the energy, also non-conserved moments, 
the stress tensor, the heat flux and fourth moments need to be in a specific form to recover 
the Navier-Stokes-Fourier equations. This was previously accomplished by assuming a pre- 
defined simple functional form, typically a polynomial, for the equilibrium population Jl| . In 
such a setting, the construction is neither unique in the choice of the discrete velocity set, nor 
in the choice of the function. Moreover, these schemes permit populations to attain negative 
values and thus make the simulation scheme unstable H |(J. The way to resolve the problem 
of non-positive form of the population is to define the equilibrium population as a minimum 
of a convex function, known as H function, under the constraint of local conservation laws. 
Recently, the advantage of such an approach was shown in the context of two dimensional 
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isothermal hydrodynamics |7], |8|, |9|, [K| . In order to avoid an explicit calculation of the local 



equilibria, these studies took an alternative root of using computationally intensive kinetic 
equations. 

Apart from the stability issue, another well known problem associated with the current 



discrete velocity models is non-adherence to the equation of the state [ XT | . In these models, 
the local equilibrium entropy does not satisfy the usual definition of the temperature as a 
function of the entropy and the energy known from the elementary thermodynamics [[Tl|1 . 

In this letter, we construct kinetic models, which are free from all the problems discussed 
above, in arbitrary dimensions. The proper choice of the set of the discrete velocities and the 
H function and the explicit expression for the equilibrium are the main result of the present 
work. These models retain the simplicity and computational efficiency of the standard 
lattice Boltzmann model. Further, for the isothermal lattice Boltzmann method an explicit 
equilibria with correct H theorem is derived. 

Before constructing such model, we briefly explain the basic setup of the discrete velocity 
models. In these methods, the kinetic equation is written for the populations /j(x, t) of 
the discrete velocities Cj, i = 1, . . . , b, defined at position x and time t. Hydrodynamic 
fields are first few moments of populations, namely p = Y^i=\ fi (density), pu a = Y^%=\ fi°ia 
(momentum density, a = 1, . . .,D, where D is the spatial dimension), and pDT + pu 2 = 
Yli=i fi c l (energy density). In the case of isothermal hydrodynamics, the hydrodynamic 
fields include p and pu a , whereas in the nonisothermal case the energy density is also included 
as independent field. Typical model kinetic equation reads, 

dth + cM^-T-'ifr-f^), (1) 

where the model collision integral on the right hand side is assumed in the Bhatnagar-Gross- 



Krook (BGK) form [12|, with r as the relaxation time. The collision integral must respect 



local conservation laws which imply (in the nonisothermal case), 

b 

^2fi*V-i ° 2 i} = {A PU a , pDT + pu 2 } (2) 

i=l 

In the isothermal case the energy constraint on the local equilibrium is excluded from this 
list. Besides Eq. (|2|), the local equilibrium must respect several other conditions for the non- 
conserved fields (examples will be given below). These latter conditions are found on the 
basis of the Chapman-Enskog analysis of the model ([!]), and they ensure that the long-time 



large-scale dynamics (|l|) would be the desired hydrodynamic equations. The problem is then 
reduced to finding a parametric expression for the equilibrium population such that all the 
constraints are satisfied. 

However, in order to construct minimal entropic model we take a different root. To begin 
with, we remind the reader of the important observation |T3|] on the relation between the 



discretization of the velocity space and the well known Grad's moment method [I4"||. Namely, 
if discrete velocities are constructed from zeros of the Hermite polynomials, the method of 
discrete velocity is essentially equivalent to Grad's moment method based on the expansion 
of the distribution function around a fixed Maxwellian distribution function. The natural 
extension of this approach towards entropic schemes is to link the discrete velocity model 
not to the Grad's method but instead to the entropic Grad's method (the maximum entropy 
approximation) [15]. To this end, Boltzmann's H function, H = J FlnFdc, where _F(x, c) 



is the one-particle distribution function, x is the position vector, and c is the continuous 
velocity, is evaluated using the Gauss-Hermite quadrature. This gives the discrete form of 
the H function, 

^,c l} = E^ ln (^)- ( 3 ) 

Here Wi is the weight associated with the ith discrete velocity Cj, and the particles mass 
and the Boltzmann constant ks are set equal to the unity. Discrete velocity populations 
/i(x) are related to values of the distribution function at the nodes of the quadrature as 
/i(x) = Wi(2 7rT )( D//2 ) exp(c^/(2T ))F(x, c,). Note that the weights are incorporated into 
the definition of fi, and that the Maxwell velocity distribution function has been factored 
out because this Gaussian probability distribution is taken into account through the Gauss- 
Hermite quadrature. Discrete- velocity entropy functions (^) for various {wi, Cj} is the unique 
input for all our constructions below. 

We shall first consider the isothermal hydrodynamics. It is known |TJ| that the mini- 
mal set of discrete velocities needed to reconstruct Navier-Stokes equations corresponds to 
zeroes of the third order Hermite polynomials in q/ \/2T . For D = 1, the three discrete 
velocities are c = {— 0, y/3 T }, whereas the corresponding weights are w = {|, |, |) 
respectively. In higher dimensions, the discrete velocities are tensor products of the discrete 
velocities in one dimension and the weights are constructed by multiplying weights associ- 
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ated with each component direction. In this way we construct the entropy function (|3]). It 
is important to remark that for D = 2 the entropy function thus obtained coincides with 
the one derived in Ref. [|TD| by a completely different kind of argument. 

The discrete-velocity local equilibrium is the minimizer of the corresponding entropy 
function under the fixed density and the momentum (Q). The explicit solution to this 
conditional minimization problem in D dimensions reads: 

3 /\/3c s 



where c s = y/Jo is the speed of the sound and the dimensionless velocity u' a = (u a /c s ). 
Note that the exponent, (<^ a /(v3cs)), in Eq. Q) takes the values ±1, and only and the 
resulting expressions for the equilibrium can be simplified in each dimensions. Equilibria 
(H) are positive definite for u a < V3c s . Without going into details of derivation, we mention 
that the equilibrium @ is the product of D one-dimensional solutions, see for example Ref. 
@. This factorization is pertinent to the derivation, and is quite similar to the familiar 
property of Maxwell's distribution function. 

We stress that none of the conditions for the higher-order equilibrium moments have 
been used while deriving (Q). Relevant higher-order moments of the equilibrium distri- 
bution, needed to establish isothermal hydrodynamics in the framework of the Chapman- 
Enskog method are the equilibrium pressure tensor, = £V /j eq Q Q Ci/3, and the equilibrium 
third-order moments, Q^L = ^ j^^ofii^c^. A direct computation shows that the present 
equilibrium results in the correct pressure tensor (correct means as obtained in the continu- 
ous kinetic theory) to the order 0(u A ), which is sufficient to simulations of the Navier-Stokes 
dynamics at small Mach number. For example, for D = 2, even at high dimensionless ve- 
locity, u a = 0.25\/T^, the error in the pressure is less than 2%. Moreover, the third-order 
moment is even more accurate as compared to the standard quadratic polynomial equilib- 
rium, which leads to a 0(u 3 ) error in simulations |16 |. In the present case, only the diagonal 



component Q^ aa have the same linear accuracy as in the standard LBM, whereas all the 
D 3 — D non-diagonal components of <5^ 7 are correct up to order 0(u 5 ). 

The set of the discrete velocities used in the present case is the same as standard D2Q9 
and D3Q27 model of the lattice Boltzmann method. The expansion of the equilibrium to the 
order 0(u 2 ) coincides with the polynomial equilibria used in the lattice Boltzmann method. 
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Thus, it is not surprising that D2Q9 model is more stable than any other LBM. However, 
due to the enhanced numerical stability and accuracy in the heat flux, the use of the positive 
definite equilibrium (Eq. ^) is more preferable in comparison to its expanded form in the 
isothermal lattice BGK method. 

In precisely the same way, the minimal entropic kinetic model for the nonisothermal 
case requires zeroes of fourth-order Hermite polynomials. For D = 1, the four discrete 
velocities and corresponding weights of Gauss-Hermite quadrature are c = (±a, ±b), and 
w = (w a ,Wb) = (T /(4a 2 ), T /(45 2 )), respectively, where a = \/3~ To) 1 / 2 , and b = 
To) 1 / 2 . The minimizer of the H function (§) corresponding to the velocity set and 
weights just described, and subject to the constraints (H), may be written as 



(fc 2 -T-^) exp{ Bc z ) if c . = ± a 

rcq I b 2 —a? exp(Ba)+cxp(-Ba) * ' 

' S P^+T-a 2 ) exp(BCi) 1fr _ + , U 

b 2 -a 2 cxp(Bb)+oxp(-Bfe) 11 Ci ^ 

Lagrange multiplier B, corresponding to the momentum constraint, has the following 
series representation: 

u 



B E 3 



3 r a 2 6 2 (a 2 + b 2 ) 



+ 0{u 5 E~ 5 ), (6) 



E A E 3 

where E = u 2 + T is the total energy density [Notice that Eq. @ is not an expansion in 
powers of velocity u, rather, in terms of u n E~ m .] Equilibrium distribution (|^) exists within 
a positivity interval, a 2 < T + u 2 < b 2 . 

In higher dimensions, the set of discrete velocities is formed by taking tensor product of 
the discrete velocities in D — 1, and the weights are products of corresponding quantities. In 
order to evaluate Lagrange multipliers in the formal solution to the minimization problem, 
/ 4 eq = iy,exp (A + B ■ Cj + Ccf), we first make an important observation that they can be 
computed exactly for u = and any temperature T within the positivity interval, a 2 < T < 
b 2 : 



b -n r - 1 ^( Wa{T ~ a2 A 

B a - U, G - TTo 775 TfrT ) 

(b 2 -a 2 ) \w b {b 2 -T)J 

Ao = l0g [ (2w a nb 2 -a 2 )° )- DaCo - 
With this, we find the equilibrium at zero average velocity and arbitrary temperature, 

6 



Ji 2 D (b 2 -a 2 ) D V) 



n 



a=l 

Factorization over spatial components is clearly seen in this solution. Once the exact solution 
for zero velocity is known, extension to u ^ is easily found by perturbation. The first few 
terms of this expansion are: 

■^'- (r-a^-r) "^ '" 4 '- 

B a = y + -^p^r {Dupueu^a^e - 3u 2 u a ) + 0(u 5 ), 
„ „ a 2 (b 2 -T)-T(b 2 -3T) 2 4 . 

For the actual numerical implementation, the equilibrium distribution function can be cal- 
culated analytically, up to any order of accuracy required, by this procedure. Note that 
dependence on the temperature in the above equations is a nonperturbative result. 

In order to establish the hydrodynamic equations corresponding to the present model, 
apart from the equilibrium pressure tensor and the equilibrium third moment <5^ T , 
one needs to check also the fourth order moment, R ^ = J2i c ia c iP c2 f^- A direct compu- 
tation shows that the equilibrium stress tensor is exact from the computational standpoint 
(accurate at least up to the order 0(u 8 )). The third moments Qa/3 7 are accurate up to the 
order 0(u6 2 ), 0(u 3 6) and 0(w 5 ), where = (To — T)/Tq is the deviation of the tempera- 
ture from the reference value. Finally, the fourth moments, R a p are accurate to the order 
O(0 2 ), O(u 2 2 ), 0(m 4 ). Thus the Navier-Stokes- Fourier dynamics is recovered to the order 
0(u 3 , 2 ). Notice that if the temperature is fixed at the reference value T = T , the moment 
<3^9 7 becomes exact to the order 0(w 5 ), unlike in the second-order accurate standard lattice 
Boltzmann models and the isothermal model constructed above. 

The local equilibrium entropy, S = —kBH^ Wt C .y(f eq ), for the nonisothermal model satis- 
fies the usual expression for the entropy of the ideal monoatomic gas to the overall order of 
approximation of the method, S = p k B In (T D I 2 / p) + 0{u 4 , 2 ). Thus, the present model is 
able to retain the thermodynamics up to the accuracy of the method. To the best of our 
knowledge, this is the first discrete velocity model which is fully consistent with thermody- 
namics. 
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When the single relaxation time BGK model ([!]) is used with the present nonisothermal 
equilibrium, the resulting transport coefficients are as follows: For D — 1, the kinematic 
viscosity v is equal to zero, while the thermal conductivity k is k = (3/2) (r p)T . For D > 1, 
we have v = (rp)T, and k = ((D + 2)/2)(rp)T [we recall that one needs to renormalize 
the relaxation time in the BGK model as r' = r p, in order to obtain density-independent 
transport coefficients]. 

Finally, we present here some details of the discretization scheme for model kinetic Eq. 
(HI). In a kinetic model there are two time scales one associated with convection and the 
other one with collisions, where the timescale of collisions is much smaller than that of 
convection. In order to have an efficient simulation scheme for the hydrodynamics, it is 
desirable to follow the time scale of the convection in the simulation. To achieve such a 
scheme we propose to simulate discrete kinetic equation 

(9) 



fi(x,5t) = L {cij5t } 



; * (X '° )+ (2T> a +8t) (f?(*°)-M*> Q )) 



were Lr Cit st} is a linear convection operator and is defined by the relation Lr c . st \ ■ g(x,t) = 
g(x — Ci5t,t), and the parameter a is defined by the condition: 

H {wuCi} (/ 4 (x, 0)) = H {wijCi} (/,(x, 0) + a (/r(x, 0) - / 4 (x, 0))) . (10) 

Close to the local equilibrium the parameter a is equals to 2St [|J. The details of the 
implementation of Eq. ( PH| ) are presented in the Ref. 0. The essence of the Eq. @ 
and ([[(]) is to separate the convection, which just shifts the population in space (operator 
L{ Ci ,st})i from the collisions which represent the relaxation of the populations towards the 
equilibrium. After each convection step, the length of collision step is restricted by the 
condition of the entropy conservation during the collision (Eq.|l0|). This lumping of many 
short collision steps ensures that rapid convergence towards hydrodynamics is achieved. The 
details of the derivation of Eq. @ and Eq. (0) from the model kinetic equation (0) will be 
presented elsewhere. 

In the isothermal case, this model can be implemented in a efficient manner by taking 



the uniform grid spacing in each direction as 5x and the time step as St = Sx/ v37o- This 
means x — Ci5t is always a grid point and the resulting method is LBM. In the nonisothermal 
case, the error in the approximations can be minimized by choosing the timestep as 5t = 
35x/b and performing a dimensional splitting, for example in the two dimensional case 

L{c lx ,c ly ,8t} ' tfOM) = L {c lx ,0,8t} " [ L {0,c ly ,5t} ' tfCM)]- 
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For directions in which velocity component equals a, we propose to discretize the one 
dimensional convection operator by the Beam- Warming operator |T7J defined as £{a,3fe/&} • 
g(x,t) = 0.02432g(x,t)+0.99784g(x-5x,t)+0.022lQg(x-25x,t). We remind the reader that 
as the CFL number approaches the value of unity (in the present case equals 3a/b = 0.95351), 
the error of the discretization vanishes [17]. The present discretization scheme, in contrast 
to the earlier proposed finite- volume and finite-difference methods |I| for solving Eq. ([I]), is 
using a large time step (of the order 0(Sx)) like LBM. 

In order to check the effectiveness of the algorithm for the nonisothermal case, we have 
performed the simulation of the Taylor vortex flow. The flow is in the isothermal and 
the incompressible set up, which is achieved in the simulation through the initial condi- 
tion. This problem is chosen to validate the theoretical expression for the viscosity and 
the discretization procedure. The flow is completely characterized by the analytical solu- 
tion u(x,y,t) = V x [(uo/k 2 ) exp [— u{k\ + k%)t] cos(fcix) cos(k 2 y)}- In the simulation, we 
have chosen uq = 0.0001, k\ = 1, k 2 = 4. The result in Fig.|Tj shows that the discretization 
procedure is working well even at very short times. 

To conclude, in this Letter we have derived minimal entropic kinetic models of the Boltz- 
mann equation for both isothermal and nonisothermal hydrodynamic simulations. The 
resulting models have correct hydrodynamics, they are equipped with the appropriate H 
function, and also they are thermodynamically consistent. A simple discretization scheme 
is proposed for the simulation. In the isothermal case, we have found analytically the cor- 
responding local equilibrium in closed form (f|), and thus proposed a isothermal LBM with 
correct H theorem. 
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FIG. 1: Simulation of the two dimension Taylor vortex flow. The velocity profiles at y = tt at three 
different times t = 0.03, 10, 50 are shown. The solid lines represent analytical results with viscosity 
v = t'T, circles represent simulation results. All the quantities are given in dimensionless units. 
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